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1 .  INTRODUCTION 


Much  progress  has  been  made  during  the  recent  years  in  the  development 
of  computational  capabilities  for  general  analysis  of  certain  nonlinear 
effects  in  solids  and  structures.  In  each  of  these  developments,  quite 
naturally,  the  first  step  was  the  demonstration  of  some  ideas  and  possi¬ 
bilities  for  the  analyses  under  consideration,  and  then  the  research  and 
development  for  reliable  and  general  techniques  was  undertaken.  This 
second  step  proved  in  many  cases  much  more  difficult,  and  in  the  case  of 
capabilities  for  analysis  of  contact  problems  has  yielded  few  general 
results. 

Although  some  of  the  first  complex  contact  problems  have  been  solved 
using  the  finite  element  method  quite  some  time  ago  [1-3],  and  much  in¬ 
terest  exists  in  the  research  and  solution  of  contact  problems  [see,  for 
ex.  refs.  4-15],  there  is  still  a  great  deal  of  effort  necessary  for  the 
development  of  a  reliable,  general,  and  cost-effective  algorithm  for  the 
practical  analysis  of  such  problems.  This  is  largely  due  to  the  fact  that 
the  analysis  of  contact  problems  is  computationally  extemely  difficult, 
despite  the  relatively  simple  mechanics  theory  used  for  these  problems. 

Much  of  the  difficulty  lies  in  that  the  boundary  conditions  of  the  bodies 
under  consideration  are  not  known  prior  to  the  analysis  but  they  depend 
on  the  solution  variables. 

The  aim  in  our  research  is  the  development  of  a  solution  algorithm 
for  analysis  of  general  contact  conditions  which  shall  include  the  possi¬ 
bilities  to  analyse, 

0  contact  between  flexible-flexible  and  rigid-flexible  bodies, 

°  sticking  or  sliding  conditions  (with  or  without  friction), 

°  large  relative  motions  between  bodies, 

°  repeated  contact  and  separation  between  the  bodies. 

Since  the  large  deformation  motion  of  the  individual  bodies  can 
in  many  cases  be  analysed  already  quite  effectively  [16],  an  algorithm 
of  the  above  nature  will  certainly  enlarge,  very  significantly,  the 
currently  available  computational  capabilities  for  practical  nonlinear 
analyses.  The  objective  in  this  paper  is  to  present  our  first  research 
results  towards  the  above  aim. 

In  this  paper  we  consider  the  large  deformation  of  two-dimensional 
planar  or  axisymmetric  bodies  in  contact  and  in  static  conditions. 

The  algorithm  we  present  contains  the  following  major  ingredients: 
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-  The  total  potential  of  the  contact  forces  is  included  in  the 
variational  formulation  to  enforce  the  geometric  compatibilities 
along  the  contact  surfaces. 

-  In  the  region  of  contact,  surface  tractions  are  evaluated  from 
the  externally  applied  forces  and  nodal  point  forces  equivalent 
(in  the  virtual  work  sense)  to  the  current  element  stresses. 

-  The  surface  tractions  between  nodal  points  (on  the  element 
segments)  are  employed  to  decide  whether  a  nodal  point  is  in 
sticking  or  sliding  contact,  or  is  releasing. 

-  The  number  of  equations  due  to  the  contact  conditions  is 
dynamically  adjusted  to  solve  two  equations  for  each  node  in 
contact  if  the  node  is  in  sticking  condition,  and  one  equation 
if  the  node  is  in  sliding  condition. 

Because  of  the  highly  nonlinear  contact  conditions  to  be  analysed, 
the  success  of  the  algorithm  largely  depends  on  an  effective  formulation 
with  special  attention  to  details.  We  believe  that  the  gradient  matrix 
used  when  sliding  conditions  are  analysed  and  the  segment  approach 
employed  to  decide  on  the  contact  conditions  are  two  important  and 
key  aspects  of  our  algorithm. 

In  the  next  two  sections  we  present  the  formulation  of  the  algorithm 
and  the  important  numerical  details.  We  have  implemented  the  solution 
method  in  the  computer  program  ADINA  [17],  and  in  Section  4  we  give 
the  solutions  to  various  sample  problems.  These  serve  to  demonstrate 
the  applicability  of  and  also  the  assumptions  used  in  the  algorithm. 


2.  FORMULATION  OF  CONTACT  PROBLEM 


Figure  1  shows  schematically  the  problem  we  consider.  This  fig¬ 
ure  shows  two  generic  bodies  which  we  arbitrarily  denote  as  contactor 
and  target.  In  the  finite  element  solution,  the  contactor  contains 
the  finite  element  boundary  nodes  that  come  into  contact  with  the 
target  segments  or  nodes.  Although  only  two  bodies  are  shown  to  come 
into  contact,  the  algorithm  can  analyse  the  contact  conditions  between 
a  number  of  bodies. 

The  basic  conditions  of  contact  along  the  contact  surfaces  are 
that  no  material  overlap  can  occur,  and  as  a  result,  contact  forces 
are  developed  that  act  along  the  region  of  contact  upon  the  target 
and  the  contactor.  These  forces  are  equal  and  opposite.  The  normal 
tractions  can  only  exert  compressive  action,  and  the  tangential  tractions 
satisfy  a  law  of  frictional  resistance. 

The  Friction  Law  Used.  Much  research  effort  is  currently  focussed 
upon  the  development  of  appropriate  friction  laws  and  the  mechanics  using 
these  laws  to  predict  motion  along  slip  surfaces,  e.q.,  [11,  18,  19]. 


PRESCRIBED 
»  FORCES 


DISPLACEMENTS 


a)  Condition  prior  to  contact 


CONTACT  REGION, 
NO  A  PRIORI 
KNOWLEDGE  OF 
REGION 


b)  Condition  at.  contact 

Figure  1.  Schematic  representation  of  problem  considered 


Considering  the  development  of  our  contact  algorithm,  we  therefore 

should  use  a  friction  model  that  is  physically  realistic  and  that 

we  can  extend  in  further  developments,  and  as  more  refined  models 

become  available.  These  criteria  are  fulfilled  using  Coulomb's  law 

of  friction  with  u  the  static  coefficient  of  friction  and  u  .  the 
s  d 

dynamic  (or  kinetic)  coefficient  of  friction  [20,  pp.  53-64]. 

Consider  the  particles  initially  in  contact:  those  belonging 
to  the  target  body  and  those  of  the  contactor.  If  tt  represents 

the  developed  tangential  tractions  along  the  contact  surfaces,  we 
assume  that  there  is  no  relative  motion  between  two  adjacent  particles 
on  the  contactor  and  the  target  in  contact,  as  long  as  | t^ |  £  u  tn, 

where  t  is  the  compressive  normal  traction  (assumed  positive).  The 

maximum  traction  of  static  friction  is  the  smallest  force  necessary 
to  start  motion.  During  motion  the  magnitude  of  the  tangential  traction 
resisted  by  friction  is  tn,  with  £  y  .  The  motion  continues 

as  long  as  the  frictional  traction  is  developed  to  equal  the  dynamic 
friction  u<j  t  that  can  actually  be  resisted.  Once  the  developed 

tangential  traction  drops  below  the  dynamic  friction,  the  relative 
motion  between  the  contactor  and  target  particles  ceases  until  such 
time  that  again  the  developed  tangential  traction  exceeds  the  fric¬ 
tional  capacity. 

We  may  note  that  with  this  friction  law,  we  neqlect  any  elasticity 
between  the  particles  in  contact  and  assume  a  rigid-plastic  contact 
behavior.  Refinements  of  this  friction  law  would  entail  the  use  of 
rate  and  state  variables,  as  discussed,  for  example  in  [19]. 

Considering  our  finite  element  formulation  of  the  above  frictional 
conditions,  we  should  note  the  following  two  important  points.  Firstly, 
although  ri oi d-pl asti c  behavior  is  assumed  between  particles  in  contact, 
the  two-dimensional  finite  element  discretization  around  the  contact 
region  can  represent  nonlinear,  e.g.,  elastic-plastic,  material  conditions. 
Secondly,  the  above  friction  law  is  in  our  finite  element  formulation 
satisfied  in  a  global  sense  over  each  individual  contact  segment  (as 
discussed  in  Section  3)  consistent  with  the  level  of  finite  element 
discretization  used. 

Some  Prel iminaries.  For  the  formulation  of  our  contact  solution 
algorithm  we  use  the  incremental  procedures  -  including  the  notation  - 
presented  in  ref.  [16,  Chapter  6]  and  recognize  that  for  each  of  the 
bodies,  the  contact  conditions  can  be  imposed  by  adding  to  the  usual 
variational  indicator,  the  total  potential  of  the  contact  forces 
with  the  constraint  of  compatible  boundary  displacements.  Hence,  in 
the  formulation  we  invoke  stationarity  of  the  following  functional. 


(1) 


7Tf  7T  -  I  W 
k  K 


where  tt  is  the  usual  (incremental)  total  potential  leading  to  the 
incremental  equilibrium  equations  without  contact  conditions,  and 

]Twk  is  the  incremental  potential  of  the  contact  forces.  This  term 

can  be  interpreted  as  a  Lagrange  multiplier  contribution  to  impose 
the  contact  conditions.  In  the  following  sections  we  concentrate 
on  the  evaluation  of  W^  for  a  generic  node  k  on  the  contactor  surface 

(and  of  the  corresponding  nodes  on  the  target  surface)  in  sticking 
and  sliding  conditions. 

Assume  that  in  the  incremental  solution,  the  response  at  time  t 
has  been  calculated  and  that  (i-1)  iterations  have  been  performed 
to  calculate  the  solution  at  time  t+At.  The  formulation  of  the 
governing  equations  is  achieved  by  establishing  for  the  next 

iteration  (i).  We  repeat  that  this  contribution  is  the  only  change 
in  the  incremental  equilibrium  equations  presented  in  ref.  [16, 
Chapter  6] . 

Figure  2  shows  a  generic  region  of  contact  considered  which 
satisfies  the  contact  conditions.  We  note  that  the  displacements 
and  coordinates  are  interpolated  linearly  between  adjacent  nodes 

on  the  contact  surfaces  of  the  bodies, ^  and  that  some  of  the 
nodes  can  be  in  contact  whereas  others  are  still  (or  again)  in 
separation.  Also,  based  on  the  assumptions  along  the  region  of 
contact,  the  contactor  nodes  cannot  be  within  the  region  of  the 
target  body,  but  the  target  nodes  can  be  inside  or  outside  the 
contactor  body.  This  point  requires  particular  attention  when 
modeling  a  problem  for  use  of  the  contact  algorithm. 

2.1  Potential  of  Contact  Forces  for  Sticking  Contact 

A  contactor  node  k  is  assumed  to  be  in  sticking  contact  under 
two  conditions: 

a)  The  contactor  node  has  penetrated  the  target  body  in  iteration 
(i-1)  whereas  it  was  not  in  contact  after  iteration  (i-2). 


^Actually,  as  will  become  apparent,  the  contact  solution  algorithm 
can  also  be  employed  when  the  contactor  and/or  target  bodies 
are  discretized  using  parabolic  elements  (see  Sections  4.1  and  4.4) 


Figure  2  Finite  element  discretization  in  contact  region. 

Nodal  point  numbers  on  contactor  surface  increase 
in  direction  such  that  when  moving  from  k  to  k+1 
the  contactor  body  is  on  the  left  hand  side. 
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b)  The  frictional  resistance  during  contact  is  sufficient  to 
prevent  sliding. 

In  case  (a)  the  contact  force  at  node  k  at  the  beginning  of  iter¬ 
ation  (i)  is  zero  and  the  contact  force  is  generated  during  iteration 
(i)  when  the  overlap  is  eliminated. 

Figure  3  shows  how  node  k  has  come  into  contact  with  the  target 
segment  j  formed  by  nodes  A  and  B,  where 


t+At  (i-1)  t+At  (i-1)  t+At  (i-1) 
-k  *  -A  ’  -B 

=  current  global  coordinates  of  nodes 
k.  A,  B,  respectively,  after  itera¬ 
tion  (i-1)  for  the  equilibrium  con¬ 
figuration  corresponding  to  time 

t+At^ 

t+At  (i-1) 

=  current  global  coordinates  of 
the  assumed  physical  point  of 
contact  of  node  k 

=  overlap 

J 

=  length  of  segment  j 

r,  s 

=  local  isoparametric  coordinate 
system  along  target  surface 

n  ,  n  =  unit  vectors  along  local  axes 
_r  -s  r,  s  on  target  segment,  respectively, 

with  respect  to  the  global  reference 
frame;  updated  during  each  iteration 
(but  for  ease  of  notation  the  super¬ 
script.  (i-1)  is  not  given) 


Note  that,  as  in  Chapter  6  of  ref.  [16],  the  left  superscript  "t+At"  on  a 
variable  denotes  the  configuration  t+At  in  the  incremental  solution,  and 
does  not  imply  a  dynamic  analysis  [16,  p.  309]. 


gure  3  Definition  of  variables  for  segment  j 


unit  vectors  along  global  x,  y 
coordinate  axes 


i,  j  ■ 


g(i-l) 


parameter  of  location  of  ohysical 
point  of  contact 


At  node  k  we  have  a  contact  force  that  we  denote  here  as 
but  whose  evaluation  we  discuss  in  detail  in  Section  3,  where 


t+At.  (i-1) 

h. 


t+At  (i-1)  _  t+At.  (i-1)  . 

-k  '  Akx  1  + 

Note  that  the  components  of 
t+At^(i-l)  Qf  Eq  (18)_ 

Also,  from  geometry. 


t+At  (i-1) 
*k 

t+At  (i-1) 

(3) 

d.n-’)  , 

J 

PrT 

(i-1)  t+At  ( 
-A 

'-'h 

(4) 

Ct+it£c( 

i-1)  t+At  (i 
-A 

-’>] 

t+At  (i-1)  . 

Aky  i 


(2) 


appear  in  the  vector 


=  R  T  rrt+At*  (i-1) 


A 


(i-1  )\  t+At  (i-1 )  -i 


And  we  have  for  target  segment  j,  the  forces  equivalent  to  - 


t+At  (i-1) 

4 


t+At,  (i-1)  _  ft( i -1 )  t+At  (i-1)  • 

ab  "  "6  h 


(7) 

(8) 


Let  the  displacement  increments  at  nodes  k,  A,  B  in  iteration  (i) 
be  Au^^,  Au^^^,  AUg^1^  respectively.  These  displacements  are  such  that 
the  overlap  A^1"^  is  eliminated.  Also,  if  contact  was  already  present, 

the  point  of  contact  C  for  node  k  is  the  same  during  each  iteration,  hence 
^  ^  =  6^  The  potential  due  to  the  contact  force  at  node  k  and  the 

corresponding  reactions  is  in  iteration  i. 


«k  *  tt4tXk(,)T  (Auk(,)  ♦  Ak(,"1))  +  tt4tXA(1)T  Aua(1) 

*tt4V1iT  AU^’  •  <9> 

Also 

tt4txk(1>  =  t+Atik(i-i)  +  tx_U)  00) 

where  AA^^  is  the  change  in  the  contact  force  at  node  k.  Using  Eqs.  (7)  to 
(10)  we  obtain 


9 


■  t*At  (1-1)T[(  (1)  t  .  (i-l);  .  4  0) 


au „(i>] 


+  AX, 


0)  [(AU  <<>  *  A  <'-'>)  -  (l-6(i-'>)  Au  .«> 


Aub'’>]  . 

This  potential  is  considered  for  all  contactor  nodes  k  that  are  in 
sticking  contact. 

2.2  Potential  of  Contact  Forces  for  Sliding  Contact 

A  contactor  node  k  is  assumed  to  be  in  sliding  contact  if  according 
to  the  criteria  given  in  Section  3,  the  tangential  force  exceeds  the 
frictional  capacity.  The  calculation  of  total  potential  for  the  sliding 
contact  condition  is  more  involved  than  for  sticking  contact  because 

the  parameter  of  location,  changes  during  iteration  (i)  to  a 

new  value  B^.  However,  the  frictional  force  is  assumed  to  remain 

constant  during  the  iteration.  Using  Eqs.  (7)  to  (9)  with  b^  we  have, 

Wk  =  t+AVl)T[(Auk(i)  +  Ak(i_1))  -  (1-B(i))  AuA(i)  -g(i)  A_uB(i)]  (12) 


whs  re 


;(i)  _  gO-1)  + 


and  from  Eq.  (5)  we  obtain,  by  1 inearization. 


4  OrT  [<AS)k(1)  *  Ak<1-'>)  -  (1-B(i-'>)ASA(1) 


AU  (’'I] 


Also,  for  sliding 

t+At,  (i)  _  t+At.  (i-1 )  .  (i) 

*k  ~  -k  A-k 

05) 

Aik(1)  =  -  AX*'1  Qs 

(16) 

where  AA^  is  the  change  in  the  magnitude  of  the  normal  component  of 

The  negative  sign  in  Eq.  (16)  is  used  because  an  increase 
in  the  normal  force  is  acting  into  the  opposite  direction  of  n^. 

Substituting  from  Eqs.  (13)  to  (16)  into  Eq.  (12)  we  obtain 


‘  ttAtXk(,',)  [<V1)+  ^0-1),  .  ayAd)  _6(i-l)  (0] 


♦  t+itx.(,-,>r(4e<1))  (aua«>  -Aub<’>) 


’ik 


+  4X^1 


-!SWU)  *  5k(i‘n)  -  O-e'1-1’)  AU„(i>] 


(17) 


where  we  neglected  (Axf 1  ))(ab^)  terms. 


This  potential  is  considered  for  all  nodes  k  that  are  in  sliding  contact. 

2.3  Governing  Finite  Element  Equations 

The  incremental  finite  element  equations  of  motion  including  contact 
conditions  are  generated  by  substituting  from  Eqs.  (11)  and  (17)  into  Eq.  (1) 
and  invoking  stationarity,  <5tt^  =  0.  Hence,  we  obtain,  using  the  usual 

procedures . 


-*  A 
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f 

t+At ^( i-1 ) 

o 

o 

g 

+ 

t+At^  (i-1) 

> 

ax^ 

< 

. 

_  ^ 

* 

-  — 

-  — 

■■  “ 

t+At? 

t+AtF (i-1) 

t+Atp  (i-1) 

- 

-c 

g 

o 

+ 

t+At.  (i-1) 

L  -c  J 

where. 


AU 


(D 


A\ 


(i) 


t+AtK(i-l ) 


t+At,,  (i-1) 
*C 


t+AtF(i-l) 


t+Atr 


t+Atp  (i-1) 

-c 


t+At  (i-1) 

-c 


NEQ 


Vector  of  incremental  displacements  in  iteration 
(i);  of  dimension  (NEQxl). 

Vector  of  increments  in  contact  forces  in  iteration 
(i);  (NEQCxl ) . 


Usual  tangent  stiffness  matrix  including  material 
and  geometric  nonlinearities  after  iteration  (i-1); 
(NEQxNEQ). 


Contact  stiffness  matrix,  for  the  effect  of  contact  conditions 
after  iteration  (i-1);  (NEQTxNEQT). 


Vector  of  nodal  point  forces  equivalent  to  element  stresses 
after  iteration  (i-1);  (NEQxl). 


Vector  of  total  applied  external  forces  at  time  t+At;  (NEQxl). 


Vector  of  updated  contact  forces  after  iteration  (i-1),  (NEQxl) 


Vector  of  overlaps  (NEQCxl). 

Total  number  of  displacement  degrees  of  freedom. 


(18) 
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NEQC  =  Total  number  of  incremental  contact  constraint  equations 

=  2x  (total  number  of  nodes  in  sticking  contact)  +  (total  number 
of  nodes  in  sliding  contact). 

NEQT  =  NEQ  +  NEQC 

Each  contactor  node  k  contributes  to  ^  and 

-c  -c 


t+.‘.t .  ( i  -1 ) 


Consider  these  terms  for  a  single  contactor  node  since  the 


contributions  for  a  number  of  nodes  are  obtained  by  addition  of  the 
individual  contributions  using  the  direct  stiffness  method  [16]. 


In  the  case  of  sticking  contact,  the  first  term  in  Eq.  (11)  results  in 
the  vector  t+AtR  whereas  the  second  term  gives  the  contact  stiffness 


-c 

t+At,,  (i-1) 


matrix  K 
-c 


,  ,  .  t+At.  (i-1) 

and  overlap  vector  A  '  ' 


t+At,  (i-1) 


t+it  (i-1) 


t+At.  (i-1) 


-d-s(f-'>)  (i-i) 

kx 

-(1-S(l_1))  t+Atx.  (i_1) 
ky 

_R(i-l)  t+At,  (i-1) 


_  (i-1)  t+At  (i-1) 


t+At.  (i-1) 


(i-1) 


V 


and  the  corresponding  solution  vector  is 


AU 


AX 


(i) 


(i) 


A“k 


A-A 


(i) 


(i) 


All, 


(i) 


L 


(i) 


AXc 


J 


Although  we  have  simply  listed  the  vectors  t+AtRc^1’^  as  shown 

in  Eqs.  (19)  and  (22),  an  important  ingredient  of  our  algorithm  is 
that  the  actual  elements  of  these  force  vectors  are  derived  as  explained 
in  the  next  section. 


Note  that  the  above  equations  (used  for  the  sample  solutions  in 
Section  4)  correspond  to  a  full  Newton  iteration.  Our  experiences  with 
the  solution  of  contact  problems  have  so  far  shown  that  for  the 
contact  equations  full  Newton-Raphson  iteration  is  usually  best. 


3.  EVALUATION  OF  STICKING  AND  SLIDING  CONDITIONS  AND  FRICTIONAL  RESISTANCE 


Much  of  the  difficulty  of  solving  contact  problems  lies  in  the 
design  of  appropriate  procedures  for  numerically  updating  the  contact 
conditions  at  a  contactor  node.  In  other  words,  the  algorithm  has  to 
decide  whether  a  node  is  not  in  contact  and  whether  the  matrices  in 
Eqs.  (19)  to  (21)  for  sticking  contact  or  the  matrices  in  Eqs.  (22)  to 
(24)  for  sliding  contact  shall  be  included  in  the  system  of  equations. 
Appropriate  decisions  during  the  iteration  concerning  these  conditions 
are  most  important  for  a  reliable  and  effective  scheme. 

After  iteration  (i-1)  the  nodal  point  displacements  t+Aty(i-l)  and 
nodal  point  forces  aR^’^  are  known  where  (see  Fig.  4) 


AR<i-l)  =  t+AtR  _  t+At  F ( i -1 ) 


(25) 
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We  note  that  at  the  nodes  not  belonging  to  a  contact  surface 

the  components  of  aR^’^  are  the  out-of-balance  loads  usually  en¬ 
countered  in  nonlinear  analysis  [16],  but  corresponding  to  the  bound¬ 
ary  nodes  affected  by  the  contact,  the  contact  forces  1 )  are 

active.  These  forces  are  evaluated  from  the  out-of-balance  aR^"^ 
and  correspond  to  tension  release,  sticking  or  sliding  conditions. 

The  procedure  of  calculating  the  contact  forces  from  the  vector 

.'.r/  1 " ^  is  effective  because  an  incrementation  of  the  Lagrange  multi¬ 
pliers  used  in  Eq.  (18)  can  -  in  other  than  geometrically  and  materially 
linear  analyses  of  sticking  contact  conditions  -  lead  to  serious  errors 
of  linearization. 

In  the  following  we  consider  the  contactor  segments  and  contactor 
nodes,  and  we  discuss  how  the  conditions  of  node  sticking  and  sliding 

can  be  reached,  and  how  the  contact  forces  are  evaluated. 

When  a  contactor  node  penetrates  the  target  body  in  an  iteration, 
which  is  decided  kinematically  by  the  displacements  of  the  contactor 
and  target  bodies  leading  to  a  geometric  overlap  (see  Fig.  3),  the 
matrices  in  Eqs.  (19)  to  (21)  are  included  in  the  solution  for  the 
next  incremental  displacements.  Hence,  in  the  first  iteration  from  no 
contact  to  a  contact  condition,  sticking  is  assumed. (+)  This  is  the 
mechanism  used  to  evaluate  the  nodal  point  forces  and  enable  a  decision 
on  whether  sticking  or  sliding  conditions  are  really  applicable. 

3 . 1  Contactor  Segment  Distributed  Tractions  and  Resultant  Forces 

The  decision  on  whether  a  contactor  node  is  releasing  or  is  in 
sticking  or  sliding  conditions  is  perhaps  most  quickly  based  on  con¬ 
sidering  the  total  and  relative  magnitudes  of  the  calculated  nodal 
points  forces.  However,  this  can  lead  to  some  numerical  difficulties, 
and  it  is  deemed  more  effective  to  establish  the  condition  at  a  con¬ 
tactor  node  from  the  accumulated  effects  and  conditions  of  the  contactor 
segments  adjacent  to  the  node. 

The  first  step  in  our  procedure  is  to  calculate  the  distribution 
of  the  tractions  along  the  contactor  boundary  given  the  nodal  point 
( i -1 )  k  k 

reactions  /.R  .  Let  t  and  t  be  the  magnitudes  of  the  distributed 

—  x  y 

tractions  (force/unit  area)  at  the  nodal  point  k  (see  Fig.  4),  then  we 
have  with  a  linear  displacement  interpolation  in  a  "consistent"  approach 
to  calculate  the  tractions  from  the  nodal  point  forces,  in  plane  stress 


We  may  note  that  for  the  special  case  of  frictionless  contact,  i.e. 
=  1.^  =  0.0,  we  can  directly  assume  perfect  sliding  conditions. 


CONTACTOR 

BODY 


a)  Out-of*balance  forces  acting  onto  contactor  body 


Figure  4.  Calculation  of  normal  and  tangential  tractions  onto 
contractor  body 


b)  Tractions  acting  onto  contactor  body 


c)  Normal  and  tangential  tractions  onto  contactor 
body.  Normal  traction  is  positive  when  acting 
onto  the  body,  tangential  traction  is  positive 
when  acting  from  node  k  to  node  (k+1). 

Figure  4  Calculation  of  normal  and  tangential  tractions  onto 
contactor  body. 


and  plane  strain  analyses,  with  uniform  thickness  h,'  ' 


where 


is  the  x-coordinate  of  nodal  point  k  at  the  end 

of  iteration  (i-1).  Using  Eqs.  (26)  and  (27)  a  tridiagonal  coefficient 
matrix  is  established  that  relates  the  unknown  tractions  to  the  known 

values  of  and  the  equations  can  be  solved  to  calculate,  in 

k  k 

each  iteration,  the  nodal  values  t  and  t  for  all  nodes  in  contact. 

x  y 

These  values  are  then  employed  to  evaluate  the  tangential  and  normal 

k  k 

segment  tractions,  t.  and  t  ,  at  the  nodes.  Note  that  to  evaluate 

1  n  k  k  k+1  k+1 

these  tractions  for  segment  j,  the  values  t  ,  t  and  t  ,  t 

x  y  x  y 

are  simply  transformed  to  the  tanqential  and  normal  directions  defined 
by  the  angle  e.  of  the  segment.  This  results  in  general  into  a  dis¬ 
continuity  of  the  normal  and  tangential  segment  tractions  at  the  nodal 
points,  see  Fig.  4. 

For  the  definition  of  the  state  of  a  segment,  we  need  the  total 
normal  and  tangential  forces  applied  to  the  seqment.  In  the  case  of 
plane  stress  and  plane  strain  analyses,  the  total  resultant  normal 

force,  TJ  ,  acting  on  seament  j  is 
n  3 

i  k+i 

V  -^(t n  *  *„  >.  <28 

and  the  total  resultant  tangential  force,  T^  ,  acting  on  segment  j  is 


,(i-l) 


V  *  V  1  • 


/  ♦  |  I 

where  d.  *  is  the  length  of  the  contactor  segment  j  in  iteration  (i), 
J 

Similarly,  for  axisymmetric  analysis  we  have 


ft+At  (i-1)  k  t+At  (i-1)  k+1, 

K  xk+l  ln  xk  tn  > 


ri 

■ 


i 

i 

i 


and 


J  = 


d<V>) 

i'2  ( 


t+Atx  (i-D  t^k  +  t+AtVi(i-l)  tM)  + 


t+At  (i-1 )  t  k  +  t+At  (i-1)  t  k+l} 
xk+l  Xk  Lt  ; 


(31) 


With  the  above  calculations  completed  the  algorithm  decides  on  the 
state  of  the  segment  using  the  segment  resultant  forces  ,  t  1  and  Coulomb' 
law  of  friction  globally  applied  over  the  segment. 

3.2  Segment  Release 

The  segment  is  assumed  to  have  experienced  tension  release  if 

is  negative,  and  in  this  case  the  contactor  segment  normal  and  tangential 
tractions  are  set  to  zero. 

3.3  Assume  Segment  was  in  Previous  Iteration  in  Sticking  Contact 

Using  the  total  normal  force  on  the  segment  TJ  the  frictional 

i  n 

capacity  of  the  segment  is  calculated  using  Coulomb's  law  of  friction, 

T^  =  u s  ,  where  us  for  the  segment  is  set  equal  to  ud  if  the  segment  was 

ever  in  sliding  (see  Section  4.2).  The  following  two  situations  can  now 
arise. 

Case  1 :  The  frictional  capacity  of  the  segment  is  larger  than  the  applied 
total  tangential  force,  i.e.,T^.  >_  |T^|  .  The  segment  continues  to  stick. 

Case  2:  The  frictional  capacity  of  the  segment  is  smaller  than  the  applied 
total  tangential  force,  i.e.,  T^  <  |T^|.  The  state  of  the  segment  is  now 

updated  to  sliding,  with  T^  =  ud  T^. 

The  results  on  whether  the  segment  continues  to  stick  or  is  now 
sliding  are  later  used  in  deciding  whether  the  contactor  nodes  are  sticking 
or  sliding  (see  Section  3.5). 

Aside  from  deciding  on  the  sticking  and  sliding  conditions  of  the 
segments,  also  the  contributions  to  the  vector  t+a'tR^('M)  must  ^ 
evaluated.  This  is  done  differently  in  the  above  two  situations. 
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In  case  1  the  distributed  tangential  and  normal  tractions  on 
the  segment  j  are  employed  to  calculate  the  nodal  point  consistent 
loads,  see  Fig.  5.  We  note  that  if  also  the  conditions  of  the  adja¬ 
cent  segments  (j-1)  and  (j+1 )  correspond  to  case  1,  these  nodal  point 
consistent  loads  are  at  the  nodes  k  and  k+1  simply  minus  the  values 


in  ',R 


(i) 


and 

t+. 


Fig.  5(b) 


summarizes  the 
Note  that  a 


In  case  2  the  segment  is  sliding, 

tractions  t^  used  in  the  calculation  of 

uniform  traction  is  assigned  such  that  in  sliding  conditions  the  total 
tangential  force  is  scaled  down  to  equal  the  frictional  capacity. 

Using  a  uniform  frictional  traction  represents  perhaps  the  simplest 
way  of  globally  satisfying  Coulomb's  law  of  friction  over  the  segment. 
Figure  5  shows  the  tractions  used  in  plane  stress  am  plain  strain 

k  +1 

analyses;  in  axisymmetric  solution  the  value  (tt  +  t^  )/2  of  Fig.  5 
is  replaced  by  t^ , 


2  V 


d(i-l)  ( t+At  ( i -1 )  t+it  (i-1) 

J  1  xk  xk+l  ) 


(32) 


In  summary,  unless  there  is  tension  release  (see  Section  3.2)  the 

contact  forces  are  calculated  as  the  consistent  nodal  point 

loads  corresponding  to  the  tangential  tractions  t^,  shown  in  Fig.  5, 

and  the  unaltered  normal  tractions  shown  in  Fig.  4.  In  tension  release 
both  the  normal  and  tangential  tractions  on  the  segment  are  set  to  zero 

3 . 4  Assume  Segment  was  in  Previous  Iteration  in  Sliding  Contact 

If  the  segment  was  in  sliding  contact,  analogous  calculations 
to  those  described  in  Section  3.3  are  performed,  but  the  friction  coef¬ 
ficient  used  is  ud.  Hence,  in  case  1  the  segment  changes  to  sticking, 

whereas  in  case  2  the  segment  continues  to  slide. 

3. 5  Conditions  of  Contactor  Nodes 

Once  the  conditions  of  the  contactor  segments  have  been  decided 

^Updating  the  tangential  tractions  in  this  way  raises  questions 
on  the  convergence  of  the  iterative  solution  as  studied  in  a 
forthcomning  communication. 


as  discussed  above,  the  algorithm  determines  the  conditions  of  the 
nodes  on  the  contactor  surface.  Table  1  summarizes  how  the  various 
conditions  (release,  sticking  and  sliding)  of  a  contactor  node  are 
reached.  We  may  note  that  these  conditions  decide  on  whether  zero 
(corresponding  to  no  contact  or  contact  release),  one  (corresponding 
to  slidinq)  or  two  (corresponding  to  sticking)  contact  equations  are 
to  be  included  in  the  incremental  equations  for  each  contactor  node. 

The  decision  on  whether  a  contactor  node  k  is  not  in  contact, 
or  is  in  slidinq  or  sticking  contact,  and  the  evaluation  of  the 
contact  matrices  (in  Eqs.  (21)  and  (24))  and  the  contact  forces  to 
be  used  in  Eq.  (18)  gives  all  the  ingredients  to  proceed  with  the 
i terat i on  ( i  ) . 


4.  SOME  SAMPLE  SOLUTIONS 


The  algorithm  presented  in  the  previous  sections  was  implemented 
in  ADI NA  and  in  the  following  i  present  the  results  of  some  sample 
analyses.  In  these  analyses,  the  primary  objective  was  to  study  the 
performance  of  the  algorithm  in  order  to  identifv  where  improvements 
are  needed  rather  than  to  solve  actual  practical  problems. 

It  is  our  experience  that  regardinq  the  solution  of  contact 
problems  some  "very  simple  looking"  problems,  including  frictional 
conditions  and  the  elasticity  of  the  structure,  may  provide  quite 
severe  tests  on  the  Performance  of  an  algorithm,  and  in  fact  may 
be  more  difficult  to  solve  than  actual  practical  engineering  problems. 

4 . 1  Analyses  of  Hertz  Contact  Problems 

Figure  6  shows  the  contact  oroblem  considered  and  the  finite 
element  idealization  used.  In  this  problem  a  long  cylinder  with 
radius  R= 1 0  was  analyzed;  hence  in  the  model  plane.strain  conditions 
wpre  assumed.  The  rigid  target  surface  was  modeled  by  specifying 
nodes  with  no  deqrees  of  freedom.  In  the  region  of  anticipated 
contact,  8-node  elements  were  used  to  model  the  contactor  with  one 
contactor  segment  always  spanning  over  one  8-node  element  side,  and 
these  elements  were  modeled  as  material lv-nonl inear-only  or  using 
the  total  Lagrangian  formulation  [16],  To  simulate  the  load  appli¬ 
cation,  the  vertical  displacements  were  prescribed  along  the  top 
' urfdce  of  the  model  and  the  total  load  P  for  a  prescribed  displace¬ 
ment  was  calculated  by  integrating  the  contact  pressures.  Figure  7 
lives  the  calculated  contact  pressures  and  a  comparison  with  the 
Hertz  analytical  solution  [21]. 

We  may  note  that  only  a  few  solution  points  are  qiven  in  Fin.  7, 
because  the  program  outputs  the  mean  traction  over  a  segment  and  at 
the  maximum  applied  load  (P  -  6600)  only  three  segments  were  in 
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TABLE  1.  STATE  OF  CONTACTOR  NODE  AS  DECIDED  BY  STATES  OF  ADJOINING  SEGMENTS 


STATE  OF  ADJOINING  SEGMENTS 

STATE  OF  NODE 

one  adjoining  segment 

other  adjoining  segment 

sticki ng 

sticking 

sliding 

tension  release 

sti  cki ng 

sliding 

s 1 i di ng 

tension  release 

sliding 

tension  release 


tension  release 


tension  release 


E  =  30000 


RIGID  TARGET  CONTACTOR 

SURFACE  SURFACE 


Ms  s  0.0 
Md  =ao 


a)  Problem  considered 

Figure  6.  Analysis  of  Hertz  plane  strain  contact  problem  (P.N.O.  and 
T.L.  formulations  denote  material ly-nonl i near-only  and  total 
Lagrangian  formulations,  resp.  [16]) 
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long  cylinder  is  modeled. 
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resp.  [16]) 
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contact.  In  order  to  obtain  more  finite  element  solution  points  and 
a  higher  solution  accuracy,  a  finer  finite  element  discretization  on 
the  contactor  surface  is  required. 

Next,  the  Hertz  contact  problem  of  a  sphere  of  radius  R= 1 00  was 
analyzed.  Hence,  Fig.  6  (a)  still  shows  the  contact  problem,  but  now 
R= 1 00 . 0  and  axisymmetric  conditions  are  considered.  Figure  8  gives 
the  finite  element  mesh  used  in  this  analysis  and, Fig.  9  shows  the 
calculated  contact  pressures  and  a  comparison  with  the  Hertz  solution. 

4.2  Motion  of  a  Rubber  Sheet  in  a  Converging  Channel 

A  sheet  of  rubber  in  plane  stress  was  confined  to  move  in  a  rigid 
horizontal  channel.  Figure  10  shows  the  sheet  and  the  finite  element 
idealization  used.(+)  The  right  face  of  the  sheet  was  subjected  to 
the  displacement  history  given  in  Fig.  11,  making  this  a  large  defor¬ 
mation  problem.  Note  that  the  displacements  were  assumed  to  vary 
slowly  so  that  inertia  effects  could  be  neglected. 

Although  the  solution  obtained  could  not  be  compared  with  an 
available  solution,  this  is  an  interesting  problem  to  study  the 
performance  of  the  contact  algorithm.  Also,  the  essential  features 
and  solution  difficulties  of  this  problem  are  frequently  encountered 
in  actual  practical  problems;  e.g.,  analysis  of  metal  forming  Drocesses. 

Figure  12  shows  the  distribution  of  normal  and  tangential  tractions 
for  different  load  steps  in  the  solution.  The  tractions  close  to  the 
face  at  which  the  displacements  are  imposed  are  not  shown  because 
a  fine  finite  element  idealization  would  be  required  to  obtain  good 
stress  predictions  near  the  face.  Note  that  the  magnitudes  of  the 
tangential  tractions,  tt>  for  times  8,  14  and  24  are  essentially 

equal  to  times  tn  —  because  practically  the  entire  rubber  sheet 

is  sliding  through  the  channel  - — and  that  the  tangential  tractions 
at  time  24  are  acting  in  the  opposite  direction  to  the  tractions  at 
times  8  and  14.  However,  at  time  18  the  tangential  tractions  have 
only  partially  reversed  and  some  segments  are  still  in  sticking  con¬ 
ditions.  It  is  this  change  in  tangential  tractions,  resulting  from 
the  reversal  in  motion,  that  is  quite  difficult  to  analyze. 

Figure  12  also  shows  the  results  obtained  when  assuming  a 
frictionless  motion.  As  expected,  for  the  frictionless  case  the 
normal  tractions  are  significantly  larger  at  times  8  and  14  (the 
imposed  displacements  increase)  and  smaller  at  times  18  and  24 
(the  imposed  displacements  decrease),  when  compared  with  the  results 
including  friction. 


^Note  that  the  rigid  target  surface  was  modeled  using  12  segments  for 
the  sole  purpose  of  demonstrating  that  contactor  nodes  can  slide 
over  target  nodes. 
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Figure  9  Solution  to  the  axisymmetric  Hertz  problem 
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a)  At  times  8  and  14 

Figure  12  Distributed  tractions  in  analysis  of  rubber  sheet  moving 
through  rigid  channel.  (Solid  line  refers  to  the 
solution  for  and  the  solution  for  u$> 

using  our  usual  algorithm;  dashed  line  refers  to  "experiment 

when  u  >  y. ) 


tn  TIME  18 


b)  At  times  18  and  24 


Fiaure  12  continued 


Finally,  the  motion  of  the  rubber  sheet  for  u  >  u.  was  also 

s  d 

analyzed  (ps  =  0.20,  ud  =  0.15).  For  this  case,  two  different  solu¬ 
tion  algorithms  were  employed.  In  the  first  solution  our  usual 
algorithm  was  used,  in  which  the  value  of  u£  for  a  segment  was  set 

equal  to  ud  as  soon  as,  and  for  all  times  thereafter,  the  static 

frictional  resistance  was  exceeded  for  the  segment.  Since  the  effect 
of  u$  >  ud  is  then  a  transient  phenomenon,  the  solution  results  for 

the  times  considered  in  Fig.  12  are  (very  closely)  the  same  as  the 
results  obtained  for  the  case  u$  =  u^.  In  the  second  analysis,  as 

an  experiment,  the  value  of  u  =  0.20  was  kept  throughout  the  solu¬ 
tion  and  the  results  marked  us  =  0.20,  =  0.15  in  Fig.  12  were 

obtained.  It  should  be  noted  that  in  this  analysis  the  time  step 

(it  =  0.5)  is  an  important  (physical)  variable  of  the  problem  — 

because  the  solution  procedure  simulates  the  frictional  motion  of 

the  sheet  for  each  solution  step  separately  —  and  it  is  questionable 

whether  the  assumptions  used  in  this  numerical  solution  appropriately 

simulate  the  actual  physical  process  of  motion.  However,  it  is 

interesting  to  note  that  for  the  chosen  values  of  u  and  relatively 

s  d  J 

small  differences  in  the  tractions  were  calculated  when  compared  with 
the  solution  for  u$  =  wd  (except  for  the  tangential  tractions  at  time 

18  which  more  drastically  changed  sign  when  u  >  uJ.  It  should  be 

emphasized  that  these  solutions  of  the  rubber  sheet  when  l  >  ... 

s  d 

should  only  be  regarded  as  a  rather  brief  numerical  experiment, 
because  there  are  many  difficult  questions  related  to  the  physics 
and  to  our  numerical  analysis  procedure  for  this  problem  that  need 
much  further  study  [19]. 

4. 3  Analysis  of  a  Snapped  Wire  in  Continuous  Writing 

The  practical  application  of  this  problem  lies  in  the  analysis 
of  the  conditions  that  arise  when  a  wire  of  a  continous  wiring  around 
a  cyl i nder  snaps . 


Figure  13  shows  the  model  considered.  Note  that  the  snapped  wire 
is  free  of  constraints  at  x  =  0  (for  y  >_  0).  The  objective  was  to 
calculate  the  stress  distribution  in  the  continuous  wiring  (modeled 
here  as  a  continuum  at  the  maximum  load  application,  i  =  300  MPa 

yy 

and  ’  =  750  MPa. 

xx 

Figure  13  (b)  shows  the  finite  element  idealization  used  for 
the  analysis  in  Fig.  14  gives  some  calculated  stresses  and  a  comparison 
with  the  results  obtained  by  Boman  [22],  The  solution  with  our  contact 


Stresses  just 


Figure  14,  continued 


algorithm  was  obtained  using  one  load  step  to  apply  the  initial 
stress  a  and  then  two  load  steps  with  At=1.0  as  shown  in  Fig.  13  (a) 

to  reach  the  final  stress  condition. 

4.4  Analysis  of  a  Buried  Pipe 

Frictional  conditions  must  frequently  be  modeled  in  the  analysis 
of  soi 1 -structure  interactions  [15,  23]. 

Figure  15  shows  a  pipe  buried  in  soil  subjected  to  the  overburden 
pressure  Pq  =  100  kPa.  The  objective  of  the  analysis  was  to  predict 

the  tractions  along  the  pipe-soil  interface.  In  this  analysis,  both 
the  pipe  and  the  soil  were  considered  linear  elastic  media,  although 
in  practice  the  soil  may  need  to  be  considered  nonlinear. 

Figure  16  shows  the  finite  element  idealization  used  for  the 
analysis,  and, fig.  17  gives  the  predicted  tractions  along  the  con¬ 
tact  surface.''  '  Also  shown  in  Fig.  17  are  the  tractions  along  the 
interface  for  the  friction  coefficients  p  =  0.0  and  p  =  °°,cal culated 
without  the  use  of  the  contact  algorithm.  These  results  have  been 
obtained  by  simply  using  constraint  equations  so  that  the  pipe  and 
soil  nodal  displacements  are  the  same  perpendicular  to  the  pipe  sur¬ 
face,  and  free  tangentially  when  p  =  0.0  and  the  same  tangentially 
when  p  =  °°. 

Figure  17  shows  that  the  results  using  the  contact  algorithm  are 
slightly  different  from  those  obtained  without  the  contact  algorithm. 
These  differences  arise  mainly  because  of  the  traction  recovery  used 
in  the  contact  algorithm,  and  because  the  contact  algorithm  always 
uses  the  updated  nodal  point  positions  (including  the  displacements) 
for  the  contact  force  calculations,  whereas  these  deformations  were 
neglected  when  using  the  constraint  equations. 


5.  CONCLUDING  REMARKS 


An  algorithm  for  the  solution  of  two-dimensional  contact  problems, 
including  large  deformation  and  frictional  conditions,  has  been 
presented.  The  solution  procedure  uses  a  Lagrange  multiplier  tech¬ 
nique  to  incrementally  impose  the  deformation  constraints  along  the 
contact  surfaces.  The  contact  forces  are  evaluated  from  distributed 
tractions  that  act  on  the  contactors.  The  tractions  are  calculated 
from  the  nodal  point  forces  (which  correspond  to  the  internal  element 
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As  seen  in  the  figure,  the  figure  p  =  »  can  be  modeled  with  the 
contact  algorithm  using  any  large  value  of  p. 
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Figure  15  Pipe  buried  in  soil  subjected  to  overburden 
pressure  Pn  =  100  kPa,  u  =  u  =  u.. 
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Figure  16  Finite  element  idealization  used  for  analysis 
of  buried  pipe 
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INTERFACE  TRACTIONS  PER  UNIT  OVERBURDEN  PRESSURE 
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Figure  17  Predicted  tractions  on  pipe/soil  interface;  solution  obtained  using 
four  load  increments  of  eoual  size 
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stresses  and  the  externally  applied  loading)  and  the  frictional 
conditions  based  on  Coulomb's  law.  The  solution  results  obtained 
using  the  algorithm  in  some  contact  problems  have  been  presented 
to  demonstrate  some  of  the  features  of  the  solution  procedure. 

Considering  the  way  frictional  conditions  are  accounted  for  in 
the  algorithm  some  important  assumptions  are  used.  First,  the 
frictional  forces  in  sliding  are  assumed  to  act  in  the  same  directions 
as  the  contact  tangential  forces  prior  to  sliding.  This  assumption 
may  require  relatively  small  load  increments  in  the  solution. 

Second,  the  frictional  calculations  make  only  use  of  the  total  tan¬ 
gential  and  normal  forces  acting  on  the  segments,  and  do  not  account 
for  the  variation  of  the  tractions  over  the  segments.  And  third, 
the  relatively  simple  friction  law  of  Coulomb  is  employed.  A  more 
refined  friction  law  would  include  rate  and  state-dependent  factors. 
Using  Coulomb's  friction  law  with  =  ud  already  a  large  number  of 

contact  problems  can  be  modeled  using  our  algorithm,  but  various 
questions  relating  to  the  physics  of  motion,  and  to  our  numerical 
solution  procedure,  must  still  be  addressed  when  us  >  wd  [19]. 

In  this  first  paper,  we  have  concentrated  on  presenting  the 
theory  used  and  on  indicating  some  applications.  Our  experiences 
with  the  algorithm  have  been  most  encouraging,  but  the  field  of 
analysis  of  contact  problems  is  very  large,  and  many  most  interesting 
aspects  relating  to  our  algorithm  deserve  further  studies,  such  as 

-  the  effect  of  the  finite  element  mesh  used  for  a  problem 
on  the  performance  of  the  solution  procedure; 

-  effective  modeling  of  the  target  and  contactor  bodies,  with 
respect  to  selection  of  an  appropriate  number  of  contactor 
and  target  segments  for  determination  of  contact; 

-  the  use  of  appropriate  load  incrementation  for  solution  of 
specific  problems,  and  the  effect  of  different  sequences 
of  load  application,  in  particular  when  ps  >  wd; 

-  the  choice  of  iteration  procedure  and  convergence  criteria 
(for  example,  perhaps  more  effective  methods  than  full  Newton 
iteration  can  be  identified); 

-  the  use  of  a  lumped  approach  for  the  traction  recovery  instead 
of  the  consistent  approach,  and  the  use  of  higher-order  contact 
segments ; 

-  rigorous  mathematical  analyses  and  convergence  studies  of 
the  algorithm  for  frictional  conditions. 

These  studies  would  be  very  valuable  because  they  will  yield 
further  insight  into  the  solution  procedure  and  provide  the  basis 
for  improvements  of  the  solution  method.  We  are  currently  pursuing 
such  studies  and  plan  to  report  upon  them  in  future  communication. 
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This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the 
reports  it  publishes.  Vour  comments/answers  to  the  items/questions  below  will 
aid  us  in  our  efforts. 
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BRL  Report  Number 
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Does  this  report  satisfy  a  need? 

(Comment  on  purpose,  related  project,  or 
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etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  _  • 

reports?  (Indicate  changes  to  organization,  technical  content,  format,  etc.) 
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